function gradient = gauss3DAffineGradient(params, mu1, sigma1, mu2, sigma2)

% compute negative gradient of closed form L2 measure between two Gaussian

dim = length(mu1);
affine = params(4:6);
trans = params(1:3);
rot = quat2rotmat(expmap2quat(expmap));
   
mu = [rot*mu1']' + trans - mu2;
sigma = rot*sigma1*rot' + sigma2;

mu_grad = gaussMeanGradient(mu, sigma, zeros(size(mu)));
sigma_grad = gaussCovarianceGradient(mu, sigma, zeros(size(mu)));

mu_rot_grad = scaleGradient(rot, mu1);
sigma_rot_grad = quadGradient(rot, sigma1);

rot_expmap_grad = rotationMatrixExponentialMapGradient(expmap);
               
rot_grad = zeros(size(rot));
for i = 1:size(rot, 1)
    for j = 1:size(rot, 2)
        m_r_g = mu_rot_grad(:, i, j);
        s_r_g = sigma_rot_grad(:, :, i, j);
        rot_grad(i, j) = sigma_grad(:)'*s_r_g(:) + mu_grad(:)'*m_r_g(:);
    end
end

trans_grad = mu_grad;
expmap_grad = zeros(3, 1);
for i = 1:3
    r_e_g = rot_expmap_grad(:, :, i);
    expmap_grad(i) = rot_grad(:)'*r_e_g(:);
end

gradient = [trans_grad; expmap_grad];
